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ABSTRACT 

We report on simultaneous broadband observations of the TeV-emitting blazar Markarian 501 be- 
t’ween 1 April and 10 August 2013, including the first detailed characterization of the synchrotron 
peak with Swift and NuSTAR. During the campaign, the nearby BL Lac object was observed in both 
a quiescent and an elevated state. The broadband campaign includes observations with NuSTAR, 
MAGIC, VERITAS, the Fermi Large Area Telescope (LAT), Swift X-ray Telescope and UV Optical 
Telescope, various ground-based optical instruments, including the GASP-WEB'T program, as well 
as radio observations by OVRO, Metsahovi and the F-Gamma consortium. Some of the MAGIC 
observations were affected by a sand layer from the Saharan desert, and had to be corrected using 
event-by-event corrections derived with a LIDAR (Light Detection And Ranging) facility. This is the 
first time that LIDAR information is used to produce a physics result with Cherenkov Telescope data 
taken during adverse atmospheric conditions, and hence sets a precedent for the current and future 
ground-based gamma-ray instruments. The NuSTAR instrument provides unprecedented sensitivity 
in hard X-rays, showing the source to display a spectral energy distribution between 3 and 79 keV 
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consistent with a log-parabolic spectrum and hard X-ray variability on hour timescales. None (of 
the four extended NuSTAR observations) shows evidence of the onset of inverse-Compton emission 
at hard X-ray energies. We apply a single-zone equilibrium synchrotron self-Compton model to five 
simultaneous broadband spectral energy distributions. We find that the synchrotron self-Compton 
model can reproduce the observed broadband states through a decrease in the magnetic field strength 
coinciding with an increase in the luminosity and hardness of the relativistic leptons responsible for 
the high-energy emission. 

Subject headings: galaxies: BL Lacs — galaxies: individual(Markarian 501) — X-rays 
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1. INTRODUCTION 

MarkarianSOl (Mrk 501) is a nearby, bright X-ray 
emitting blazar at z = 0.034, also known to emit very- 
high-energy (VHE ; E > 100 GeV) gamma-ray photons 
(|Quinn et al.111995) . Blazars are among the most extreme 
astrophysical sources, displaying highly variable emission 
at nearly every wavelength and timescale probed thus 
far. These objects are understood to be active galactic 
nuclei that are powered by accretion onto supermassive 

Kavli Institute for Cosmological Physics, University of 
Chicago, Chicago, IL 60637, USA 

School of Physics and Center for Relativistic Astrophysics, 
Georgia Institute of Technology, 837 State Street NW, Atlanta, 
GA 30332-0430 

Department of Physics and Astronomy, Barnard College, 
Columbia University, NY 10027, USA 

Department of Physics and Astronomy, University of 
California, Los Angeles, CA 90095, USA 

Physics Department, California Polytechnic State Univer¬ 
sity, San Luis Obispo, CA 94307, USA 

Department of Applied Science, Cork Institute of Technol¬ 
ogy, Bishopstown, Cork, Ireland 

Argonne National Laboratory, 9700 S. Cass Avenue, 
Argonne, IL 60439, USA 

Max-Planck-Institut fiir Radioastronomie, Auf dem Huegel 
69, 53121 Bonn, Germany 

Institut de Radio Astronomie Millimetrique, Avenida 
Divina Pastora 7, Local 20, 18012 Granada, Spain 

Institute of Astronomy, Bulgarian Academy of Sciences, 72 
Tsarigradsko shosse Blvd., 1784 Sofia, Bulgaria 

Centre for Space Research, Private Bag X6001, North-West 
University, Potchefstroom Campus, Potchefstroom, 2520, South 
Africa 

Graduate Institute of Astronomy, National Central Uni¬ 
versity, 300 Zhongda Road, Zhongli 32001, Taiwan 

Astronomical Observatory, Volgina 7, 11060 Belgrade, 
Serbia 

Istanbul University, Science Faculty, Department of As¬ 
tronomy and Space Sciences, Beyazi t, 34119, Istanbul, Turkey 
Aalto University Metsahovi Radio Observatory, 
Metsahovintie 114, FI-02540 Kylmala, Finland 

Institute of Astronomy and NAO, Bulgarian Academy of 
Sciences, 1784 Sofia, Bulgaria 

Department of Physics, Brigham Young University Provo, 
UT 

School of Cosmic Physics, Dublin Institute For Advanced 
Studies, Ireland 

Institute for Astrophysical Research, Boston University, 
725 Commonwealth Avenue, Boston, MA 02215 

Astronomical Institute, St. Petersburg State University, 
Universitetskij Pr. 28, Petrodvorets,198504 St. Petersburg, 
Russia 

Research Institute for Science and Engineering, Waseda 
University, 3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan 

Abastumani Observatory, Mt. Kanobili, 0301 Abastumani, 
Georgia 

Engelhard! Astronomical Observatory, Kazan Federal 
University, Tatarstan, Russia 

Aalto University Metsahovi Radio Observatory, 
Metsahovintie 114, 02540 Kylmala, Finland 

Aalto University Department of Radio Science and Engi¬ 
neering, P.O. BOX 13000, FI-00076 AALTO, Finland 

Institute of Astronomy with NAO, BAS, BG-1784, Sofia, 
Bulgaria 

^^Astron. Inst., St.-Petersburg State Univ., Russia 
Pulkovo Observatory, St.-Petersburg, Russia 
Isaac Newton Institute of Chile, St.-Petersburg Branch 
INAF, Osservatorio Astronomico di Torino, 10025 Pino 
Torinese (TO), Italy 

Department of Physics, University of Colorado Denver 
Denver, CO 

Astronomical Observatory, Volgina 7, 11060 Belgrade, 
Serbia 

European Southern Observatory, Karl-Schwarzschild-Str. 
2, 85748 Garching, Germany 


black holes and hav e relativistic jets pointe d along the 
Earth’s line of sight (lUrrv &: Padovanill 199511 . Relativis¬ 
tic charged particles within blazar jets are responsible 
for the non-thermal spectral energy distribution (SED) 
which is characterized by two broad peaks in the vF^ 
spectral representation. The origin of the lower-energy 
peak is relatively well understood, resulting from the syn¬ 
chrotron radiation of relati vistic leptons in the presence 
of a tangled magnetic field (lMarschei]l2008ll . Within the 
leptonic paradigm, the higher-energy SED peak is at¬ 
tributed to inverse-Compton up-scattering by the rela¬ 
tivistic leptons within the jet of either the synchrotron 
photons themse lves, namely synchro tron self-Compton 
(SSC) emission (iMaraschi et al.l 119921) . or a photon field 
external t o the jet, namely ext e rnal Compton (EC ) emis¬ 
sion (e.g. IDeriner et al.l flQ^ ISikora et al.l 1199411 . Al¬ 
ternatively, hadronic models attribute the higher-energy 
peak of blazar emission to proton synchrotron emis¬ 
sion and/or synchrotron emissi on by secondary lep¬ 
tons produced i n p-y interactions (|Aharonian et al.l[2002l : 
lBednareEll993ll . 

Along with the other nearby VHE blazar Mrk 421, 
Mrk 501 represents one of the most comprehensively 
studied VHE blazars. The blazar has been the sub- 
ject of multiple broa d band observation campaigns (e.g. 
Cataiiese et al.l Il997t iKataoka et all 119991 : iPetrv etahl 
20nnHAbdo et al.ll2nilall . Mrk 501 is one of the bright¬ 
est X-ray sources in the sky, and has been observed by 
RX TE to display signif icant X-ray variability up to 20 
keV (IGliozzi et al.ll200^ . During a phase of high activity 
at VHE energies in 1997, this source was also observed 
by BeppoSAX to display unusually hard, correlated X- 
ray e mission up to >1 00 keV, with a photon index of 
r < 2 (IPian et al.lll'^ . 

Observations of Mrk 501 have so far lacked sufficient 
sensitivity at the hard X-ray energies (10-100 keV). Ob¬ 
servations at hard X-ray energies provide direct insight 
into the highest energy particles through detection of 
synchrotron emission. There is also the possibility for 
insight into the lower energy particles through the de¬ 
tection of inverse-Compton emission from photon up- 
scattering by the lower-energy electrons. As a relativistic 
synchrotron emitter, the falling edge of the synchrotron 
peak mimics the energy distribution of the emitting par¬ 
ticles, allowing the highest energy particles to be directly 
probed through hard X-ray observations. The energy- 
dependent cooling timescale can lead to more rapid vari- 
abilit y at hard X-ray en ergies than at soft X-ray ener¬ 
gies. Iciiozzi et al.l (|2006[1 reported independent soft (2- 
10 keV) and hard (10-20 keV) X-ray variability of Mrk 
501 using RXTE. 

Other hard X-ray observat ions have previously been 
performed wit h BeppoSAX (IMassaro et aLll^OOdall and 
Suzaku HXD (|Anderhub et al.ll2009li . Due to the rapid 
X-ray variability displayed by blazars such as Mrk 501, 
the long integration time required for significant detec¬ 
tion and spectral reconstruction by the aforementioned 
X-ray instruments was not ideal for extracting informa¬ 
tion about hard X-ray variability. Much more sensitive 
hard X-ray observations of blazars, however, are now pos¬ 
sible with Nuclear Spectroscopic Telescope Array NuS- 
TAR. 

NuSTAR is a hard X-ray (3-79 keV) observa¬ 
tory launched into a low Earth orbit in June 2012 
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(|Harrison et al.l[2013f) . It features the first focusing hard 
X-ray telescope in orbit that allows high sensitivity be¬ 
yond the 10 keV cutoff shared by all other currently 
active focusing soft X-ray telescopes. The inherently 
low background associated with concentrating the X-ray 
light enables NuSTAR to achieve approximately a one- 
hundred-fold improvement in sensitivity over the colli¬ 
mated and coded-mask instruments that operate in the 
same spectral range. 

NuSTAR observed MrkSOl four times in 2013 as part 
of a simultaneous inultiwavelength (MWL) campaign, in¬ 
cluding VHE observations by MAGIC and VERITAS, 
high-energy (HE; 100 MeV-100 GeV) gamma-ray obser¬ 
vations by the Fermi Large Area Telescope (LAT), soft 
X-ray and UV observations with Swift X-ray Telescope 
(XRT) and Ultraviolet Optical Telescope (UVOT), opti¬ 
cal observations from a number of ground-based instru¬ 
ments including the GASP-WEBT program, as well as 
radio observations by the Owens Valley Radio Obser¬ 
vatory (OVRO; 15 GHz), Metsahovi (37 GHz) and the 
F-Gamma monitoring program, providing measurements 
between 2.64 GHz and 228.39 GHz. The NuSTAR obser¬ 
vations took place on 2013 April 13, 2013 May 8, 2013 
July 12 and 13 (MJD 56395, 56420, 56485 and 56486, 
respectively), with the latter two observations resulting 
from target of opportunity (ToO) exposures triggered by 
an elevated state observed by the Swift XRT and the 
MAGIC telescopes. 

We use these observations to study the hard X-ray 
spectral behavior of Mrk 501 in detail over multiple 
flux states. The NuSTAR observations, analysis and re¬ 
sults are detailed in Section 2, with the contemporaneous 
MWL observations, analysis and results shared in Sec¬ 
tion 3. After comparing the simultaneous Swift XRT and 
NuSTAR observations in Section 4, we investigate vari¬ 
ability of the source in Section 5. The MWL SEDs are 
constructed over the multiple observed states and inves¬ 
tigated in terms of a single-zone equilibrium synchrotron 
self-Compton model in Section 6, with discussion and 
conclusions provided in Section 7. 

2. NUSTAR OBSERVATIONS AND ANALYSIS 

In order to maximize the strictly simultaneous overlap 
of observations by NuSTAR and ground-based VHE ob¬ 
servatories during this broadband campaign of Mrk 501, 
the observations were arranged according to visibility of 
the blazar at the MAGIC and VERITAS sites. The NuS¬ 
TAR coordinated observations involving both VERITAS 
and MAGIC were performed on 2013 April 13 and 2013 
May 8, with the NuSTAR ToO observations (initiated by 
Swift and MAGIC) performed on 2013 July 12 and 13. 
The NuSTAR observations typically spanned 10 hours, 
resulting in 10-30 ks of source exposure after removing 
periods of orbital non-visibility. The observation details 
are summarized in Table 1. The data were reduced using 
the standard NuSTARDAS software packag(0 vl. 3.1. 

The spectral analysis was performed with XSPEC 
Version 12.7.1. The data were binned to require 20 
counts per bin, and fit with three spectral models via 
minimization. The first model applied to the data is 
a power law 






Figure 1. The spectral energy distributions of Mrk 501 derived 
from the Nu-STAR observations, showing the PL (red), BKNPL 
(green) and LP (blue) models fitted to each observation. The NuS¬ 
TAR observations show significant detection of the blazar up to at 
least 65 keV in each observation. The data-to-model ratios are 
shown in the bottom panel of each plot, with the spectral fit pa¬ 
rameters summarized in Table 2. Spectra have been rebinned for 
figure clarity. 

A(E)pl = KiE/Eo)-^, (1) 

referred to as the PL model for the remainder of this 
work, where E{E) is the flux at energy E, F is the index, 
K is the normalization parameter (in units of photons 
“^s“^) and Eq is fixed at 10 keV. 

applied to the data is a 


'http://heasarc.gsfc.nasa.gov/docs/nustar/aiialysis/ keV ^Cin 

' http://heasarc.nasa.gov/docs/software/lheasoft/xanadu/xspec/X^]jgM^^^^j^gd;| model 
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Table 1 

Summary of the NuSTAR hard X-ray observations of Mrk 501. The observations are sometimes referred to with the last three digits of 

the Observation ID within this work. 


Observation 

ID 

MJD Exposure 

Range 

Exposure 

[ks] 

Number 

Orbits 

Detection 
Range [keV] 

60002024002 

56395.1-56395.5 

19.7 

6 

3-60 

60002024004 

56420.8-56421.5 

28.3 

10 

3-65 

60002024006 

56485.9-56486.2 

11.9 

4 

3-70 

60002024008 

56486.8-56487.1 

11.4 

4 

3-70 


Table 2 

NuSTAR spectral fit summary, with integral flux values (in units of xl0“^^ erg cm~^ derived from the log-parabolic fits. Data, 

models and ratios are shown in Figure 1. The indices of the LP fits are derived at 10 keV. The errors for each parameter are found using 
a value of Ax^=2.706, corresponding to a 90% confidence level for a parameter. Observation IDs are shortened by removing the first 

60002024 identifier in column one. 


Obs. 

ID 

Power law 


Broken Power law 



Log Parabola 


Index 

r 

PL 

X^/DOF 

Index 

Fi 

Index 

V2 

Tbreak 

[keV] 

BKNPL 

xVdof 

Index 

r 

Curvature 

d 

LP 

X^/DOF 

3-7 keV 
Flux 

7-30 keV 
Flux 

002 

2.216±0.009 

831/700 

2.04±0.03 

2.34±0.02 

6.3±0.4 

747/698 

2.290±0.010 

0.26±0.03 

729/699 

3.72±0.02 

4.81±0.03 

004 

2.191±0.006 

1204/889 

1.25±0.20 

2.21±0.01 

3.1±0.1 

1211/887 

2.250±0.008 

0.21±0.02 

1051/888 

5.19±0.02 

6.98±0.05 

006 

2.060±0.006 

1246/924 

1.92±0.02 

2.22±0.02 

7.9±0.4 

1057/922 

2.115±0.008 

0.24±0.02 

1024/923 

12.08±0.09 

18.6±0.1 

008 

2.081±0.007 

1152/863 

1.90±0.02 

2.25±0.02 

7.4±0.3 

914/861 

2.149±0.008 

0.32±0.02 

892/862 

10.75±0.05 

16.4±0.1 


broken power law, referred to as BKNPL model for the 
remainder of this work. The model is made up of two 
power-law photon indices, meeting at a break energy 

k 

^(^^)bKNPL = (2) 

where Fi and r 2 represent the photon indices below and 
above the break energy Flbreak, respectively. 

The third spectral model applied to the data is a log 
parabola, referred to as the LP model for the remainder 
of this work. This model has been suggested to bet¬ 
ter represent the X-ray spectra o f TeV-detected blazars 
between 0.2 and 100 keV (e.g. IMassaro et al.l I2004bl : 
iTramacere et al.|l2007all . This model allows the spectral 
index to vary as a function of energy according to the 
expression 

A{E)i^P = (3) 

with a curvature parameter (3. The spectral data, model 
fits and data-to-model ratios for each NuSTAR observa¬ 
tion are shown in Figure 1. The spectral fitting results 
for each model as applied to the NuSTAR observations 
are summarized in Table 2. The errors for each parame¬ 
ter are found using a value of Ax^=2.706, corresponding 
to a 90% confidence level for one parameter. 

For all four NuSTAR observations, the X-ray emission 
of Mrk 501 is bes t represented with a lo g parabola. A sta¬ 
tistical F—test (|Snedecor et al.l [198311 using the and 
degrees of freedom (DOF) of the PL versus LP fit re¬ 
sults in F-statistics of 97.8, 129.3, 200.1 and 251.3 for 
the observations 002, 004, 006 and 008, respectively, 
corresponding to probabilities of 1.1x10“^^, 4.6x10“^®, 
2.9x10“^^ and 7.9x10“®° for being consistent with the 
null PL hypothesis. The broken power-law fit to the 
second NuSTAR observation, ID 004, produces a break 
energy at the lower limit of the NuSTAR sensitivity win¬ 
dow, and is interpreted as a failed fit. The other three 


observations fit the break energy near Fbreak=7 keV, mo¬ 
tivating the decision to present the NuSTAR flux values 
in the 3-7 keV and 7-30 keV bands throughout this work. 
The upper bound of 30 keV is the typical orbit-timescale 
detection limit for the Mrk 501 observations. 

The NuSTAR observations show the blazar to be in 
a relatively low state for the first two observations, and 
a relatively high state during the last two observations, 
with the 3-7 keV integral fluxes derived from the log- 
parabolic fits 2-4 times higher than found for the first 
two observations. More specifically, the average 3-7 
keV integral flux values (in units of 10“^^ erg cm“^ 
s“^) were 3.72±0.02 and 5.19±0.02, respectively, for the 
observations occurring on MJD 56395 and 56420, and 
12.08±0.09 and 10.75±0.05, respectively, for the observa¬ 
tions starting on MJD 56485 and 56486. In the same flux 
units, the 7-30 keV integral flux values for the first two 
observations are similarly 3-4 times lower than the flux 
states observed in the last two observations (4.81±0.03 
and 6.98±0.05 on MJD 56395 and 56420 as compared to 
18.6±0.1 and 16.4±0.1 on MJD 56485 and 56486). These 
integral flux values are summarized in Table 2. 

The NuSTAR observations extend across multiple oc- 
cultations by the Earth, and the integral flux and index 
(F) light curves for the orbits of each extended observa¬ 
tion are shown in Figure 2. The periods with simultane¬ 
ous observations with the ground-based TeV instruments 
of MAGIC and VERITAS are highlighted by grey and 
brown bands in the upper portion of each light curve. 
The observations and results from MAGIC and VERI¬ 
TAS for these time periods are summarized in Section 
3.1. 

The 3-7 keV and 7-30 keV integral flux values of the 
first exposure (Observation ID 002) show low variability 
(x^ = 7.0 and 13.4 for 5 DOF), while the trend of in¬ 
creasing flux in both the 3-7 keV and 7-30 keV bands 
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Figure 2. The NuSTAR orbit-binned light curves, with 3-7 keV 
(black) and 7-30 keV (grey) integral flux values (top panel of each 
plot) and the log-parabolic indices (F, lower panel) with the curva¬ 
ture parameters (/3) fixed to the value found for the full NuSTAR 
exposure. The third and fourth observations are shown in the third 
plot. The periods where simultaneous quality-selected observations 
with MAGIC and VERITAS occurred are highlighted in the top 
panel of each plot with color coded bands. We note that the ver¬ 
tical axes are set differently for each observation to allow a clear 
view of the orbit-to-orbit variability and that the light curve for 
the full campaign is shown in Figure 5. 

is clear during the second observation (Observation ID 
004). The 7-30 keV flux increases from (5.1±0.1) xl0“^^ 
erg cm“^ s“^ to (8.8±0.1) xl0“^^ erg cm“^ in fewer 
than 16 hours. The 7-30 keV increases from (1.7±0.1) 
xl0“^° erg cm“^ s“^ to (2.0±0.1) xl0“^° erg cm“^ s“^ 
in fewer than 7 hours on MJD 56485 (Observation ID 
006) and significantly decreases from (1.9±0.1) xl0“^° 
erg cm“^ to (1.4±0.1) xl0“^° erg cm“^ s“^, again 
in fewer than 7 hours on MJD 56486 (Observation ID 
008). 

The relation between the log-parabolic photon indices 
and 7-30 keV flux values resulting from the fits to the 
NuSTAR observations of Mrk 501 are shown for each 
observation separately in Figure 3. The curvature /3 was 
not seen to change significantly from orbit to orbit and 



Flux [ ergs cm'^ s'^ ] 


Figure 3. The log-parabolic fit index F at 10 keV versus the 
7-30 keV integral flux for NuSTAR, binned by orbit. The first 
exposure is shown in red, the second in violet, and the last two in 
cyan, with solid lines meant to guide the eye along the parameter 
evolution over the full observations. In all three cases, the spectrum 
hardens when the intensity increases; in the fourth observation, the 
spectrum then softens as the intensity decreases. 

therefore was fixed at the average value found for each 
observation (see Table 2 for values). The count rate light 
curves show no indications of variability on a timescale 
of less than an orbit period (^^90 minutes ). As observed 
previ ously in the X-ray band for Mrk 501 (|Kataoka et al.l 
Il999f) . the source was displaying a harder-when-brighter 
trend during this campai gn. This has also been observed 
in the past for Mrk 421 (jT^kahashi et al.lll996l) . 

3. BROADBAND OBSERVATIONS 
3.1. Very-High-Energy Gamma Rays 
3.1.1. MAGIC 

MAGIC is a VHE instrument composed of two imaging 
atmospheric Cherenkov telescopes (lACTs) with mirror 
diameters of 17 m, located at 2200 m above sea level at 
the Roque de Los Muchachos Observatory on La Palma, 
Canary Islands, Spain. The energy threshold of the sys¬ 
tem is 50 GeV and it reaches an integral sensitivity of 
0.66% of the Crab Nebula flux above 220 GeV with a 
50-hour observation (jAleksic et al.ll^lSaT) . 

MAGIC observed Mrk 501 in 2013 from April 9 (MJD 
56391) to August 10 (MJD 56514). On July 11 (MJD 
56484), ToO observations were triggered by the high 
count rate of ^15 counts s“^ observed by Swift XRT (see 
Section 3.3). The flaring state was observed intensively 
for five consecutive nights until July 15 (MJD 56488). 
After that the observations continued with a lower ca¬ 
dence until August 10. 

The source was observed during 17 nights, collecting a 
total of 22 hours of data with zenith angles between 10° 
and 60°. Only five hours survived the standard quality 
cuts for regular MAGIC data analysis because many ob¬ 
servations were taken during the presence of a Saharan 
sand-dust layer in the atmosphere known as “Calima”. 
As we explain below, using the LIDAR information we 
could recover 10 of the 17 hours which would have been 
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Figure 4. MAGIC and VERITAS spectra averaged over epochs 
with simultaneous NuSTAR exposures. The power-law spectral 
fitting parameters for the VHE data are summarized in Table 3. 
Only statistical (Icr) error bars are shown for each of the spectral 
points. 

rejected otherwise. The telescopes were ope rated in the 
so-called wobble mode (|Fomin et alJ I1994I1 , where the 
pointing direction is changed every 20 (or 15) minutes 
among 2 (or 4) positions with an offset of 0.4° from the 
source position. 

All th e data were analyze d following the standard pro¬ 
cedure (|Aleksic et al.ll2C)l^ using th e MAGIC Analysi s 
and Reconstruction Software lMAR.S: |Zanin et al.ll2013l l. 
An image cleaning was applied based on information 
of signal amplitude and timing of each pixel, and the 
shower im ages were p arametrized using the Hillas pa¬ 
rameters (|Hillad [1985D . For the reconstruction of the 
gamma-ray direction and the gamma-hadron separation, 
the random forest method is applied using t he image pa¬ 
rameters and the stereos copic parameters. (lAlbert et al.l 
120081 : lAleksic et nil201fl( l. The energy reconstruction uti¬ 
lizes look-up tables. The analysis steps were confirmed 
independently with data from the Crab Nebula and ded¬ 
icated Monte Carlo simulations of gamma-ray showers. 

A fraction of the dataset (10.4 of 15.1 hours, specif¬ 
ically the observations between MJD 56485 and MJD 
56514) was affected by “Calima,” a Saharan sand-dust 
layer in the atmosphere. A correction within the frame¬ 
work of the MARS software is applied to account for 
the absorption due to Calima using LIDAR measure¬ 
ments taken_smmlta^ously with the MAGIC observa¬ 
tions (jFruck et al.ll201,‘tl . The correction was carried out 
in two steps. Due to the dust attenuation during Cal¬ 
ima, the estimated energy is shifted towards low ener¬ 
gies, and thus is corrected event by event, as the first 
step. Then, to account for the shift of the energy esti¬ 
mation, a correction to the collection area is applied as a 
second step, due to the energy dependence in the collec¬ 
tion area. The atmospheric transmission values for this 
method were obtained from the temporally closest LI¬ 
DAR measurement. During the observations affected by 
Calima the atmospheric transmission ranged from 85% 
down to 60%, being relatively stable within a timescale 


of one day, which is a typical feature of a Calima layer 
(unlike a cloudy sky). The precision on the energy cor¬ 
rection is estimated to be around 5% of the attenuation 
(40% to 15%), which corresponds to < 2% of the esti¬ 
mated energy, at most. After the Calima correction, the 
energy threshold increases inversely proportional to the 
transmission value. This correction method was tested 
independently on a Crab Nebula d ataset observed unde r 
similarly hazy weather conditions (iFruck fc Gaua 1201511 . 
Details of the method can be found jn jFruckl (j2014 ). This 
is the first time an event-by-event atmospheric correction 
is applied to MAGIG data. 

The analysis results of the MAGIC data taken dur¬ 
ing good weather conditions have a systematic uncer¬ 
tainty in the flux normalization and in the energy scale. 
For both of them, the component changing run-by-run 
is est imated to he ^11% using Crab Nebula observa¬ 
tions (jAleksic et al.l[^I5all . It is attributed mainly to the 
atmospheric transmission of the Cherenkov light, which 
can change on a daily basis (even during so-called good 
weather conditions) and the mirror reflectivity, which can 
change also on a daily basis due to the deposition of dust. 
The atmospheric correction applied in the analysis of the 
data taken during Calima increases this run-by-run sys¬ 
tematic error from 11% to 15% due to the uncertainty in 
the correction. Since the systematic uncertainty can be 
different according to the atmospheric correction,we have 
added 15% or 11% (with or without the atmospheric cor¬ 
rection) to the statistical errors of the flux in quadrature 
for the evaluation of flux variability. 

The summary of the MACIC analysis results for obser¬ 
vations occurring simultaneously with NuSTAR is pro¬ 
vided in Table 3. The derived spectra are shown in Fig¬ 
ure 4, where the spectral points are drawn with statis¬ 
tical errors only. The resultant flux values above 200 
GeV range from (2.39 ±0.51) x 10“^^ ph cm“^ s“^ (0.11 
Grab Nebula flux) on MJD 56395 to (5.52±0.87) x 10-i° 
ph cm“^ s“^ (2.5 times the Grab Nebula flux) on MJD 
56484. As seen in the overall light curve (top panel of Fig. 
5, shown again only with statistical errors), MAGIC ob¬ 
servations indicate a significant variability around MJD 
56484. A hint of intra-night variability was observed on 
MJD 56486 and 56487 simultaneously with the NuSTAR 
observations, as shown in the zoomed-in light curve (top 
panel of Figure 6). During these two nights the VHE 
emission is consistent with a constant flux, resulting in a 
X^/DOF of 7.3/4 (12% probability) with the inclusion of 
the systematic error. Without accounting for the addi¬ 
tional systematic error, the constant fit to the flux results 
in a xVDOF of 57/4. 

3 . 1 . 2 . VERITAS 

VERITAS is a VHE instrument comprised of four 12- 
m lACTs and is se nsitive to gamma rays between ~100 
GeV and -30 TeV (iHolder et al.l[200l:IKiedall2013ll . This 
instrument can detect 1% Crab Nebula flux in under 25 
hours. VERITAS observed Mrk 501 fourteen times be¬ 
tween 2013 April 7 (MJD 56389) and 2013 June 18 (MJD 
56461), with 2.5 and 1.0 hours quality-selected exposures 
occurring simultaneously with NuSTAR on MJD 56395 
and MJD 56421, respectively. On days without simulta¬ 
neous NuSTAR observations, the exposure times ranged 
between 0.5 hours and 1.5 hours. The observations oc- 
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Table 3 

MAGIC and VERITAS observations, analysis and spectral fit summary for AuSTAiJ-simultaneous observations. Observations occurring 
on the same day are grouped with horizontal lines. Daily average values of MAGIC observations are shown in bold, below the results for 
each observation occurring on that day. Statistical (Icr) error bars are provided for the power-law indices and the integral fluxes. The flux 
value between MJD 56486.106 and 56486.148 (shown in italics) is estimated with fitting parameters due to an energy thresh old above 200 
GeV. The significance of the observed gamma-ray signals is computed according to Eqn. 17 in lLi &: Mai II1983I1 . 


Exposure 
Start MJD 

Exposure 
Stop MJD 

Exposure 

Length 

[hr] 

Instrument 

Zenith 

angle 

[deg] 

Detection 

Significance 

Power-law 

Index 

Integral Flux 
> 200 GeV 
[xl0“^^ ph cm“^s“^] 


DOF 

56395.179 

56395.223 

1.0 

MAGIG 

10-14 

7.8 

2.50±0.24 

2.39±0.44 

0.58 

6 

56395.336 

56395.493 

2.5 

VERITAS 

15-35 

8.3 

3.1±0.4 

1.85±0.38 

0.76 

5 

56421.142 

56421.209 

1.1 

MAGIG 

12-28 

12.5 

2.24±0.08 

5.08±0.54 

15.5 

13 

56421.340 

56421.462 

1.0 

VERITAS 

20-32 

14.7 

2.25±0.15 

4.45±0.61 

6.9 

9 

56485.972 

56486.014 

1.0 

MAGIG 

12-24 

20.4 

2.19±0.07 

20.8±1.2 

10.0 

12 

56486.039 

56486.083 

1.0 

MAGIG 

28-43 

20.7 

2.39±0.08 

25.2±1.3 

26.5 

10 

56486.106 

56486.148 

1.0 

MAGIG 

48-60 

14.3 

2.71±0.12 

32.4±2.0 

11.9 

11 

56485.972 

56486.148 

2.9 

MAGIC 

12-60 

32.3 

2.28±0.04 

24.3±0.8 

24.1 

15 

56486.966 

56487.022 

1.3 

MAGIG 

12-27 

25.2 

2.37±0.06 

24.9±1.1 

20.3 

12 

56487.050 

56487.091 

0.9 

MAGIC 

33-46 

18.5 

2.23±0.09 

17.8±1.0 

14.5 

11 

56486.966 

56487.091 

2.2 

MAGIC 

12-46 

31.8 

2.31±0.05 

20.9±0.7 

30.4 

12 


curring simultaneously with NuSTAR are summarized in 
Table 3. Due to an annual, ~2 month long monsoon sea¬ 
son in southern Arizona where VERITAS is located, no 
VERITAS observations were possible for this campaign 
after 2013 June 18. 

The VERITAS observations were taken with 0.5° off¬ 
set in each of the four cardinal dire ctions to enable s i- 
multaneous background estimation (jFomin et al.l 11993) . 
Events were reconstructed fo llowing the procedure out¬ 
lined in lAcciari et al.l (j2008a[) . The recorded shower im¬ 
ages were parameterized by their principal moments, giv¬ 
ing an efficient suppression of the far more abundant 
cosmic-ray background. Cuts were applied to the mean 
scaled width, mean scaled length, apparent altitude of 
the maximum Cherenkov emission (shower maximum), 
and 0, the angular distance between the position of Mrk 
501 and the reconstructed origin of the event. The results 
were independently reproduced w ith two analysis pack¬ 
ages (iCoganI [2008t lPrc)koDbll2013D . The uncertainty on 
the energy calibration of VERITAS is estimated at 20%. 
Additionally, the systematic uncertainty on the spectral 
index is estimated at 0.2, app earing to be relat ively in¬ 
dependent of the source slope (|Madhavanll2013ll . 

A differential power law is fit to the data {dN/dE oc 
E~^) to characterize the VHE spectrum of the source. 
VERITAS observed Mrk 501 to vary by no more than a 
factor of three in flux throughout the observations, with 
the integral flux ranging from (1.85±0.38) x 10“^^ ph 
cm“^s“^ above 200 GeV (8% Crab Nebula flux above the 
same threshold) on MJD 56395 to (4.45±0.61) x 10“^^ 
ph cm“^s“^ (20% Crab Nebula flux) on MJD 56421. The 
source displayed low spectral variability, ranging between 
F = 3.1 ± 0.4 in the low flux state to F = 2.19 ± 0.07 
in the higher flux state. The observation and analysis 
results are summarized in Table 3 (for NuSTAR simul¬ 
taneous observations only), with the VHE spectra of the 
NuSTAR simultaneous observations shown in Figure 4. 
Day-to-day uncertainties in flux calculations that might 
be introduced by different atmospheric conditions (even 
under strictly good weather conditions) are not included 
in Table 3 and are estimated at less than 10%. 

3.1.3. VHE Results 


The full light curve of VHE observations from MAGIC 
and VERITAS is shown in Eigure 5, with a zoom into 
the period of elevated flux in Figure 6. The flux val¬ 
ues are shown with statistical errors only. The MAGIC 
and VERITAS observations of Mrk 501 in 2013 show the 
source in states which are consistent with the range of 
states observed in the past. The observations of VERI¬ 
TAS, occurring primarily in the beginning of the cam¬ 
paign, detected the source in a 5-10% Crab state, in 
agreement with the early MAGIC observations. Later 
on in the campaign, MAGIC observed a flux elevated 
state of order ~ 2.5 times the Grab flux. 

3.2. High-Energy Gamma Rays 

Fermi LAT is a pair-conversion telescope sensitive 
to photons between 30 MeV and several hundred GeV 
(I Atwood et al.l[20?)^ . Spectral analysis was completed 
for two periods contemporaneous with the NuSTAR 
observations using the unbinned maximum-likelihood 
method implemented in the LAT ScienceTools soft¬ 
ware package version v9r31pl, which is available from 
the Eermi Science Support Center. The LAT data be¬ 
tween MJD 56381 and MJD 56424 was used for compar¬ 
ison with the first two NuSTAR exposures, while MJD 
56471 to MJD 56499 was used for NuSTAR exposures 
occurring during the elevated state. 

“Source” class events with energies above 100 MeV 
within a 12° radius of Mrk 501 with zenith angles < 100° 
and detected while the spacecraft was at a < 52° rock¬ 
ing angle were used for this analysis. All sources within 
the region of interest from th e second Fermi LAT cata¬ 
log f2FGL. lNolan et aI~ll2012D are included in the model. 
With indices held fixed, the normalizations of the compo¬ 
nents were allowed to vary freely during the spectral fit¬ 
ting, which was performed using the instrument response 
functions P7REP_S0URCE_V15. The Galactic diffuse emis¬ 
sion and an isotropic component, which is the sum of the 
extragalactic diffuse gamma-ray emission and the resid¬ 
ual charged particle background, were modeled using the 
recommended files 0 The flux values were computed us- 

The files used were gll_ieni_v05_revl. f it for 

the Galactic diffuse and iso_source_v05.txt for 
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Figure 5. The broadband light curves of Mrk 501 from MJD 56380 to 56520. The VHE data are shown with statistical error bars only. 
Optical data are corrected as described in Section 3.4. All radio light curve points for 2-llOmm are provided by the F-Gamma consortium. 


ing an unbinned maximum likelihood analysis while hx- 
ing the spectral indices for the sources within the region 
of interest. The systematic uncertainty of the LAT ef¬ 
fective area is estimated as 10% below 100 MeV and de¬ 
creasing linearly in Log(E) to 5% between 316 MeV and 

10 GeVra 

The light curve for LAT observations of Mrk 501 was 
computed between MJD 56380 and 56520 in week-long 
bins (second panel from the top in Figure 5) and 3.5-day 
bins between MJD 56474 and 56488 (second panel from 
top of Figure 6). Single day-binned light curve was also 
investigated, but no day within the time period provided 


a signiHcant detection. More specihcally, n o day pro¬ 
vided a test statistic iTSi lMattox et al1ll996ll of greater 
than 9. 

During the Hrst epoch (MJD 56381-56424), the spec¬ 
tral analysis of the LAT data shows the blazar had an in¬ 
tegral flux of Fo,i-iooGeV=(5.3±4.4)xlO“®ph cm“^s“^, 
and an index ofT = 2.0±0.3. Analysis of the 
second epoch (MJD 56471-56499) results in an inte¬ 
gral flux of Fo,i-iooGeV=(6.5±2.1)xlO“®ph cm“^s“^ 
and index ofT = 1.7 ±0.1. These values are con¬ 
sistent with the average flux and index values calcu¬ 
lated over the first 24 months of the science phase 
of the LAT mission and reported in the 2FGL cata- 

“^s“^ and r = 


the isotropic diffuse component, both available at 
http: //f ermi. gsf c. nasa. gov/ssc/data/access/lat/BackgroundModels 
http://fermi.gsfc.nasa.gov/ssc/data/analysis/ 

LAT.caveats.html 


(Fn.i-in nnpv=(4.8ibl . 9)xl 0 cm 

±0.03: lN^n et al. llBTl . 
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The broadband light curve zoomed in to the period of the elevated X-ray and VHE gamma-ray state. 


Figure 6. 

3.3. Swift X-Ray and UV Telescope Observations 

The XRT onboard Swift (iGehrels et al.l[200^ is a fo¬ 
cusing X-ray telescope sensitive to photons with energies 
between 0.3 and 10 keV. The Swift satellite observed Mrk 
501 59 times between 2013 January 1 and 2013 Septem¬ 
ber 5 (MJD 56293 to 56540). All XRT observations were 
carried out using the Windowed Timing (WT) readout 
mode. The data set was first processed with the XRT- 
DAS software package (v.2.9.0) developed at the ASI Sci¬ 
ence Data Center and distributed by HEASARC within 
the HEASoft package (v. 6.13). Event files were cali¬ 
brated and cleaned with standard filtering criteria with 
the xrtpipeline task using the calibration files as available 
in the Swift CALDB version 20140120. 

The spectrum from each observation was extracted 
from the summed and cleaned event file. Events for 


the spectral analysis were selected within a circle of 20 
pixel 46") radius, which encloses about 80% of the 
Swift XRT point spread function (PSF), centered on the 
source position. The background was extracted from a 
nearby circular region of 40 pixel radius. The ancillary 
response files were generated with the xrtmkarf task, ap¬ 
plying corrections for PSF losses and CCD defects using 
the cumulative exposure map. The latest response ma¬ 
trices (v.014) available in the Swift CALDB were used. 
Before the spectral fitting, the 0.3-10 keV source energy 
spectra were binned to ensure a minimum of 20 counts 
per bin. 

The data were fit with an absorbed power-law model, 
with index F, as well as an absorbed log-parabolic model, 
where in both cases the neutral hydrogen column density 
was se t at 1.55 xlO^^cm”^, taken from iKalberla et al.l 
(|2005f) . The summary of the XRT observations and spec- 
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3-7 keV Swift XRT-Measured Flux [ erg cm'^ s‘^ ] 


Figure 7. The power-law index versus 3-7 keV flux values fit to 
the Swift XRT observations of Mrk 501. 

tral analysis results are provided in Table 4. The light 
curve of the observations, including 0.3-3 keV and 3-7 
keV integral flux bands, is shown in Figure 5, with a 
zoom into the period of elevated flux in Figure 6. The 
3-7 keV band is not traditionally quoted for Swift XRT 
data, but is motivated by direct comparison to the 3-7 
keV band computed for the NuSTAR observations. 

Mrk 501 displays a relatively steady flux state until 
after MJD 56480, when the flux increases to (38.3±1.5) 
xl0“^^ ergs cm“^s“^ on MJD 56483 (corresponding to 
the day with the XRT count rate of 15 counts s“^ which 
triggered MAGIC and NuSTAR observations). This high 
X-ray state was followed by a general drop in flux, contin¬ 
uing through the last XRT observation included in this 
work (2013 September 1; MJD 56540). 

The power-law fitted indices and 3-7 keV flux de¬ 
rived from the power-law fits are plotted in Figure 7 
for all 59 observations. The source clearly displays the 
harder-when-brighter trend fou nd previously in othe r 
TeV blazars, such as Mrk 421 (|Takahashi et al.lll996f) . 
This behavior is similar to that displayed in the hard 
X-ray band 7-30 keV observed by NuSTAR and shown 
in Figure 3. Notably, the photon indices in the soft X- 
ray band are systematically harder than those observed 
by NuSTAR in the 7-30 keV band. The spectral index 
observed by Swift XRT (F, determined at 1 keV) ranges 
between 1.4 and 2.2 (Figure 7) while the NuSTAR index, 
determined at 10 keV, ranges from 2.1 to 2.4 (Figure 3). 

Additionally, UV/optical observations were collected 
with the UVOT onboard Swift. These observations were 
carried out using the “filter of the day”, i.e. one of the six 
lenticular filters (V, B, U, UVWl, UVM2, and UVW2), 
unless otherwise specified in the ToO request, so images 
are not always available for all filters. There are 50 obser¬ 
vations included in this Mrk 501 campaign, 18 of which 
included exposures in all filters while the remaining 32 
observations contain UV imaging only. 

For each filter observation, we performed aperture pho¬ 
tometry analysis using the standard UVOT software dis¬ 


tributed within the HEAsoft 6.10.0 package and the cali¬ 
bration included in the latest release of CALDB. Counts 
were extracted from apertures of 5" radius for all fil¬ 
ters and co nverted to f l uxes using the standard zero 
points from iPoole et al.l (j2008[l . The flux values were 
then de-reddened using the value of E{B — V) = 0.017 
(iSchaflv fc Finkbeinerl I 20 T 1 II with A\/E{B — V) ratios 
calculated for UVOT filters usin g the mean Galact ic in¬ 
terstellar extinction curve from iFitzpatricld (|1999ll . No 
variability was detected to occur within single exposures 
in any filter. The processing results were verified, check¬ 
ing for possible contamination from nearby objects falling 
within the background apertures. 

3.4. Optical 

Temporal coverage at optical frequencies was pro¬ 
vided by various telescopes around the world, includ¬ 
ing t he GASP-WEBT program le.g. iVillata et al.ll2008L 
l2009f) . In particular, we report observations performed in 
the i?-band from the following observatories: Crimean, 
Roque de los Muchachos (KVA), Lulin (SET), Abas- 
tumani (70cm), Skinakas, Rozhen (60cm), Vidojevica 
(60cm), Perkins, Liverpool, St. Petersburg, West Moun¬ 
tain Observatory (WMO), the robotic telescope net¬ 
work AAVSOnet, the 60 cm and 1 m telescopes at 
the TUBITAK National Observatory (TUG TOO and 
TUG TlOO) and the Fred Lawrence Whipple Observa¬ 
tory (FLWO). Host galaxy es t imati on for the R filter 
is obtained from lNilsson et all (j2007l l. with apertures of 
7.5" and 5" , used for the various instruments. Galactic 
extin ction was accounted for acco rding to the coefficients 
from iSchaflv fc F inkbein ^ (120111) . The calibration stars 
reported in iVillata et al.l(jl998| l were used for calibration. 

Due to different filter spectral responses and anal¬ 
ysis procedures of the various optical data sets (e.g. 
for signal and background extraction) in combination 
with the strong host galaxy contribution (^^12 mJy 
for an aperture of 7.5" in the i?-band), the reported 
fluxes required instrument-specific offsets of a few mJy. 
These offsets are introduced in order to align multi¬ 
instrumental light curves, and were determined using 
several of the GASP-WEBT instruments as reference, 
and scaling the other instruments using simultaneous ob¬ 
servations. The required offsets for each instrument are 
as follows: Abastumani (70cm)=4.8 mJy; Skinakas=1.2 
mJy; Rozhen (60cm)=-1.3 mJy; Vidojevica (60cm)=2.2 
mJy; St.Petersburg=0.3 mJy; Perkins=0.6 mJy; Liver- 
pool=0.6 mJy; AAVSOnet=-3.4 mJy; WMO= -0.7 mJy; 
TUG T60=0.5 mJy; TUG T100=-1.2 mJy. Addition¬ 
ally, a point-wise fluctuation of 0.2 mJy (~0.01mag) was 
added in quadrature to the statistical errors in order to 
account for potential differences of day-to-day observa¬ 
tions within single instruments. Within Figure 5, the 
i?-band observations can be seen to remain fairly steady 
around 4.5 mJy. 

3.5. Radio 
3.5.1. Metsdhovi 

The 14-m Metsahovi Radio Observatory also partici¬ 
pated in this multi-instrument campaign, as it has been 
doing since 2008. Metsahovi observed Mrk 501 every few 
days at 37 GHz. Details of th e observing strategy an d 
data reduction can be found at iTerasranta et al.l (Il998f) . 
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Table 4 

Swift XRT observations and analysis results for A^iiSTjdiJ-simultaneous periods. Integral flux values are calculated according to the PL 
model, and are provided in xlO~^^ erg cm~^ s“^ units. The errors for each parameter are found using a value of Ax^=2.706, 

corresponding to a 90% confidence level for a parameter. 


Observation 

ID 

Date 

[MJD] 

Exp 

[s] 

Flux 
2-10 keV 

Flux 

0.5-2 keV 

Flux 

3-7 keV 

Flux 

0.3-3 keV 

Index 

r 

xVdof 

r 

LP 

0 

LP 

X^/DOF 

00080176001 

56395.06 

9636.0 

6.9±0.1 

6.41±0.06 

3.6±0.1 

ll.OiO.l 

2.05±0.01 

403.5/416 

2.06±0.02 

-0.02±0.04 

402.6/415 

00091745001 

56485.84 

250.7 

21.1±1.7 

12.7±0.4 

10.9±0.9 

22.3±0.7 

1.77±0.05 

108.1/94 

1.74±0.08 

0.10±0.16 

107.0/93 

00030793235 

56485.98 

709.1 

24.3±1.1 

14.6±0.2 

13.1±0.9 

24.1±0.4 

1.77±0.03 

228.7/222 

1.75±0.05 

0.03±0.09 

227.6/221 

00030793236 

56486.31 

1002.0 

24.0±0.7 

14.1±0.3 

13.4±0.6 

23.4±0.4 

1.73±0.03 

291.6/270 

1.68±0.04 

0.13±0.08 

285.1/269 

00030793237 

56487.04 

949.5 

19.1±0.9 

12.0±0.2 

10.4±0.4 

18.9±0.3 

1.76±0.03 

229.9/237 

1.73±0.05 

0.07±0.08 

228.9/236 


As can be seen in the bottom panel of Figure 5, there 
is evidence of a low level of variability at 37 GHz as ob¬ 
served by Metsahovi. This variability is quantified in 
terms of fractional variability (see Section 5.1). 

3.5.2. OVRO 

Regular 15 GHz observations of Mrk501 were carried 
out using the O VRO 40-m telescope w ith a nominal bi¬ 
weekly cadence (jRichards et al.l[2011[ l. The instrument 
consists of off-axis dual-beam optics and a cryogenic high 
electron mobility transistor low-noise amplifier with a 
15 GHz center frequency and 3 GHz bandwidth. The 
two sky beams were Dicke-switched using the off-source 
beam as a reference, while the source was alternated be¬ 
tween the two beams in an ON-ON mode to remove at¬ 
mospheric and ground contamination. The total system 
noise temperature was about 52 K. The typical noise level 
achieved in a 70-second observation was 3-4 mJy. The 
flux density uncertainty includes an additional 2% un¬ 
certainty mostly due to pointing errors, but does not 
include the systematic uncertainty in absolute calibra¬ 
tion of about 5%. Calibration was performed using a 
temperature-stable diode noise source to remove receiver 
gain drifts; the flux density s cale is derived from observa¬ 
tions of 3C 286 assuming the iBaars et'aTI 1)197711 value of 
3.44 Jy at 15 GHz. Details of t he reduction and calib ra- 
tion procedure can be found in lRichards et al.l l)2011[ l. 

3.5.3. F-Gamma 

The cm/mm radio light curves of MrkSOl were ob¬ 
tained within the framework of a Fermi-related monitor- 
ing program of g a mma- r ay blazars (F-Gamma p rogram; 
iFiihrmann et al.l (|2007[ 1: lAngelakis et al.l (|2008[ 1'1. The 
millimeter observations were closely coordinated with the 
more general flux monitoring conducted by IRAM, and 
data from both programs are included here. The overall 
frequency range spans from 2.64 GHz to 142 GHz using 
the Effelsberg 100-m and IRAM 30-m telescopes. 

The Effelsberg measurements were conducted with the 
secondary focus heterodyne receivers at 2.64, 4.85, 8.35, 
10.45, 14.60, 23.05, 32.00 and 43.00 GHz. The obser¬ 
vations were performed quasi-simultaneously with cross¬ 
scans; that is, slewing over the source position, in az¬ 
imuth and elevation direction with an adaptive number 
of sub-scans for reaching the desired sensitivity (for de¬ 
tails, see Fuhrmann et al. 2008; Angelakis et al. 2008). 
Subsequently, pointing offset correction, gain correction, 
atmospheric opacity correction and sensitivity correction 
were applied to the data. 


The IRAM 30-m observations were carried out with 
calibrated cross-scans using the Eight Mixer Receiver 
(EMIR) horizontal and vertical polarization receivers op¬ 
erating at 86.2 and 142.3 GHz. The opacity-corrected 
intensities were converted to the standard temperature 
scale and finally corrected for small remaining pointing 
offsets and systematic gain-elevation effects. The con¬ 
version to the standard flux density scale was done us¬ 
ing the instantaneous conversion factors derived from fre¬ 
quently observed primary (Mars, Uranus) and secondary 
(W3(OH), K3-50A, NGG7027) calibrators. 

4. SIMULTANEOUS NUSTAR AND SWIFT EXPOSURES 

Since Mrk 501 is highly variable, detailed inferences re¬ 
garding the broadband SED and its temporal evolution 
require simultaneous observations of multiple bands. In 
particular, for the determination of the low-energy peak 
Fsyn, and the flux at Fgyn, F(Fsyn), Swift XRT and NuS- 
TAR observations must be simultaneous. There are five 
periods within the campaign for Mrk 501 where the ob¬ 
servations by NuSTAR and Swift occurred within one 
hour of each other. The Swift exposure IDs for these 
quasi-simultaneous periods are summarized in Table 4. 
For Mrk 501, Egyn is located in the X-ray band and can 
be determined reliably (except for the first NuSTAR ob¬ 
servation where Egyn is < 0.85 keV) since there is no 
evidence of X-ray variability of Mrk 501 on a time scale 
shorter than a NuSTAR orbit (~ 90 minutes). 

As a precursor to the joint fitting of XRT and NuS¬ 
TAR data, we confirm agreement between the 3-7 keV 
flux values derived from the Swift XRT and NuSTAR 
fitted models. There is a residual discrepancy (not a 
uniform offset) at the level of < 10%. Using XSPEC, we 
performed simultaneous fitting to the datasets using the 
absorbed log-parabolic model as done in Section 2 for 
the NuSTAR data alone. During the fitting process, we 
allowed the normalizations of the data sets to vary, but 
required the same spectral shape parameters. A repre¬ 
sentative plot of the simultaneous fit for XRT and NuS¬ 
TAR data collected on MJD 56485 is provided in Figure 
8. The model spectrum is shown as a solid line in Fig¬ 
ure 8. The agreement between XRT and NuSTAR was 
stu died and found to be within the calibration uncertain- 
tie d^^H 

For the determination of the spectral parameters char¬ 
acterizing the synchrotron peak (namely the energy Egyn 
and F{Esyn)) with the simultaneous NuSTAR and Swift 

http://heasarc.gsfc.nasa.gov/docs/heasarc/caldb/swift 
/docs/xrt/SWIFT-XRT-CALDB-09_vl8.pdf 
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Figure 8. Example of a broadband X-ray spectrum of Mrk 501 in the crucial region where the synchrotron peak (in the E X F{E) 
representation) is located. The spectra result from a simultaneous observation with Swift (green) and NuSTAR (FPMA: red, FPMB: 
black) on 2013 July 12-13. The spectral fit used a log-parabolic model (see the text) with Galactic column density of 1.55 X 10^^ cm“^. 
For the purpose of illustrating the intrinsic spectrum of the source, the solid lines which represent the fit to the Swift and NuSTAR data 
show the spectrum before the Galactic absorption. The normalizations of the Swift and NuSTAR data were allowed to be free, and the 
offset between them was less than 10%, thus illustrating generally good cross-calibration of the two instruments. 


XRT observations, we apply the log-parabolic model 
modified by the photoelectric absorption due to our 
Galaxy, with a (fixed) neutral hydrogen column density 
of 1.55 X 10^° cm“^, taken from iKalberla et al.l (j2005D . 
The procedure to search for Egyn involves the variation 
of the “normalization energy” parameter (in the logpar 
model in XSPEC) until the local index T returns a value 
of 2 — then i?syn corresponds to the peak in the E xF{E) 
representation. This procedure correctly accounts for the 
effect of the soft X-ray absorption by Galactic column 
density as the absorption is included in the model fitted 
to the data. For the determination of the error on Esym 
we freeze the “local index” — defined at energy Egyn — 
to a value of 2, and then step the value of Egyn keeping 
all other parameters free. We then search for the value of 
the F'yjj which corresponds to the departure of from 
the minimum by = 2.7. The error quoted is the dif¬ 
ference between Egyn and The E^yn and curvature 

parameters (/3) for each of the simultaneous data sets are 
summarized in Table 5. We quote the value of F{Esyn) 
inferred from the NuSTAR module FPMA (Focal Plane 
Module A). 

The combination of Swift XRT and NuSTAR observa¬ 
tions provides an unprecedented view of the synchrotron 


peak variability. From Table 5, it is evident that the syn¬ 
chrotron peak moves by a factor of about ten during this 
campaign, with the highest synchrotron peak occurring 
during the elevated X-ray and gamma-ray state. 

5. VARIABILITY 
5.1. Fractional Variability 

In order to quantify the broadband variations we utilize 
the fract ional variability , Fyar. We follow the description 
given in IVaughan et al.l (120031) , where Avar is calculated 
as: 




^2 - (a2) 

(F,)^ 


(4) 


where (Fy) is the average photon flux, S is the standard 
deviation of the flux measurements, and (cr^) is the mean 
squared error of the measurement. 

Fvar was determined for the temporal binning and sam¬ 
pling presented in Figure 5 and Table 3 (for MJD 56485 
and 56486, the bold lines in Table 3 are used). The value 
of Fvar is known to be dependent on sampling and should 
be interpreted with caution. For example, a well sam¬ 
pled light curve with small temporal bins will allow us to 














































































































































































14 


Furniss et al. 


Table 5 

Fitting results for Swift XRT and NuSTAR simultaneous observations. The data were simultaneously fit with a log-parabolic function. 


Observation 

ID 

Date 

[MJD] 

Orbit 

Number 

-F'syn 

[keV] 

F{Eayii) 

[xl0“^^ ergs cm“^s“^] 

Curvature 

/3 

X^/DOF 

60002024002 

56395.1 

1 

<0.85 

4.1 

0.061 

669/673 

60002024006 

56485.9 

1 

4.9±0.7 

13.8 

0.21 

596/577 

60002024006 

56486.0 

2 

5.1±0.9 

13.7 

0.22 

697/715 

60002024006 

56486.2 

4 

7.0±0.8 

14.6 

0.2 

877/848 

60002024008 

56487.1 

4 

3.3±0.9 

11.2 

0.17 

832/851 
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Figure 9. The fractional variability (FVar) calculated for each 
instrument separately. 

probe the variability on small timescales (e.g. NuSTAR), 
which could be hidden if the variability is computed with 
fluxes obtained with relatively coarse temporal bins (e.g. 
Fermi LAT). 

The fractional variability for each band (from 15 GHz 
radio through VHE) is shown in Figure IHl For the pe¬ 
riod of observations covered in this work, the fractional 
variability shows a double-peaked shape with the high¬ 
est variability in the X-ray and VHE bands. A similar 
broadband v ariability pattern has recently been reported 
for Mrk501 (LDoerf |Ml3; Aleksic et al.ll2015cf). for Mrk 
421 (|Aleksic et al.ll2015l: : Balokovic et al.l 1201511 and for 
other high-synchrot ron-peaked blazars in, for example, 
lAleksic et all (|2014ll . This double-peaked shape of Evar 
from radio through VHE can be interpreted as resulting 
from a correlation between the synchrotron and inverse- 
Compton peaks. 

Fvar is below ~5% at 15 GHz and optical/UV frequen¬ 
cies, while at 37 GHz the fractional variability is ^20%. 
The relatively high fractional variability at 37 GHz is not 
produced by any single flaring event, but rather by a con¬ 
sistent flickering in the radio flux. Such flickering is not 
typically o bserved in bla z ars. bu t has been reported for 
Mrk 501 in lAleksic et al.l (|20i5dl . At X-ray frequencies, 
Evar gradually increases with energy, reaching the largest 
value (^0.6) in the 7-30 keV band measured by NuSTAR. 
The Evar computed for the Swift XRT 3-7 keV observa¬ 
tions is higher than for the NuSTAR 3-7 keV fluxes due 
to the larger temporal coverage of the Swift observations, 
allowing for observation of Mrk 501 during high activity 
levels that were not observed with NuSTAR. 

The Swift X RT Evar for Mrk 501 published in 
iStroh fc Faicond (|2013ll was 0.15 or 0.18, depending on 


the timescale used for calculation, illustrati ng that the 
value o f Evar is dependent on sampling. In lAbdo et al.l 
(|2011all . RXTE-ASM (2- 10 keV) and Swift BAT (15- 
50 keV) show Evar values between 0.2 and 0.3, although 
it should be noted that due to the limited sensitiv¬ 
ity of RXTE-ASM and Swift BAT (in comparison with 
Swift XRT and NuSTAR), the variability was studied on 
timescales larger than 30 days. 

5.2. Cross Correlations 

Cross-correlations between the different energy bands 
were studied with the Discrete Correlatio n Function 
(DCF) described in lEdelson fc KroIiS (jl988[l . The DCF 
method can be applied to unevenly sampled data, and no 
interpolation of the data points is necessary. Also, the 
errors in the individual flux measurements are naturally 
taken into account when calculating the DCF. One im¬ 
portant caveat, however, is that the resulting DCF versus 
time lag relation is not continuous, and hence the results 
should only be interpreted with a reasonable balance be¬ 
tween the time resolution and the accuracy of the DCF 
values. It is also important to only consider instruments 
with similar time coverage. In this study, we considered 
all the energy bands with a non-zero fractional variabil¬ 
ity. Among the Swift UVOT data, only the UVW2 filter 
was checked, as it is the filter which has the best time 
coverage across the Swift UVOT observations and also is 
least contaminated by the host galaxy light. For a better 
time coverage, MAGIC and VERITAS data points are 
combined to make a single data set as the VHE band. 

A significant correlation in the DCE was seen only be¬ 
tween the VHE data and the 0.3-3 keV and the 3-7 keV 
Swift XRT bands. Eor both of the combinations, the 
largest correlation is seen with a time lag of 0 ± 1.5 days. 
This result does not change if the binning of 3 days is 
altered. Note that the NuSTAR observations covered a 
relatively short period with a dense sampling, thus we did 
not see any significant correlation between NuSTAR and 
any other band. Since the observations of Swift XRT 
and NuSTAR were made simultaneously (within a few 
hours) with the VHE observations, correlations between 
the X-ray and the VHE observations were investigated 
in more detail (see Section 5.3). R 

5.3. X-ray/VHE Correlation 

The light curve of the broadband observations is shown 
in Eigure 5, with a zoom of the period showing an ele¬ 
vated X-ray and VHE state in Eigure 6. The VERITAS 
and MAGIC flux points within the light curve are shown 
with statistical errors only. Correlation studies using the 
VHE flux values are completed with statistical and sys¬ 
tematic errors included, as described below. The radio, 
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optical and UV observations show relatively steady flux 
over the campaign period, while the largest amplitude of 
variability can be seen in the X-ray and VHE gamma- 
ray bands. An elevated state in both the X-ray and VHE 
bands can be seen to occur on MJD 56483 {Swift Obser¬ 
vation ID 00030793232 in Table 4). Zooming in on this 
epoch (Figure 6 ), shows that the NuSTAR observations 
occurring on MJD 56485 and 56486 occurred after the 
highest state observed by MAGIC and Swift. The XRT 
observations show an elevated X-ray flux in both the 0.3- 
3 and 3-7 keV bands on MJD 56483. 

A comparison between the NuSTAR-ohseived X-ray 
photon flux values (derived from XSPEC) in the 3-7 and 
7-30 keV bands and the epochs of simultaneous VHE ob¬ 
servations is shown in Figure 10. During this campaign, 
10 observations occurred within one hour between either 
NuSTAR and MAGIC (7 observations) or NuSTAR and 
VERITAS (3 observations). The simultaneous X-ray and 
VHE data, where the VHE data include both statistical 
and systematic errors, were fit with both a linear and a 
quadratic function. 

Within the one-zone SSC emission paradigm, there is a 
physical motivation for a quadr atic relationship betwee n 
the X-ray and VHE flux values (jMarscher fc Geai]ll985l l. 
More specifically, the inverse-Compton flux depends not 
only on the density of photons, but also on the density 
of the electron population producing those photons. If, 
however, the particle population is energetic enough for 
the inverse-Compton scattering to occur in the Klein- 
Nishina regime, the relationship between the X-ray and 
VHE fluxes can be complex and will depend in detail 
on the energy bands considered, the particle energy loss 
mecha nisms and the mag netic field evolution. In par¬ 
ticular, [Katarzyhskl (j2005h suggest that a roughly linear 
relationship may arise during the declining part of a flare 
when the emitting region expands adiabatically, leading 
to a decrease of both the particle number density and 
the magnetic field strength. 

A quadratic relationship provides a better fit than the 
linear fit for the 3-7 keV flux values measured simultane¬ 
ously by NuSTAR, with of 11.4 and 87.3, respectively, 
for 9 DOF. The 3-7 keV flux and the > 200 GeV flux 
are highly correlated, with a Pearson correlation coeffi¬ 
cient (r) of 0.974. Similarly, for the 7-30 keV band, the 
quadratic relation fits the data better than the linear re¬ 
lation, with of 17.5 and 79.1, respectively, for 9 DOF. 
The r-value for the 7-30 keV flux and the > 200 GeV 
flux is 0.979. 

A comparison between the Swift-ohsevved X-ray pho¬ 
ton flux values (derived from XSPEC) in the 0.3-3 and 
3-7 keV bands and the epochs of simultaneous VHE ob¬ 
servations is shown in Figure 11. These data are not 
simultaneous with the NuSTAR observations shown in 
Figure 10 and therefore the results cannot be directly 
compared. During this campaign, 12 absolutely simulta¬ 
neous observations occurred between Swift and MAGIC 
(10) and Swift and VERITAS (2), shown in Figure 11. 
Similarly as done for the NuSTAR bands, the simulta¬ 
neous Swift X-ray and VHE data were fit with both a 
linear and a quadratic function with an offset fixed to 
zero. For the C).3-3 keV flux values measured simultane¬ 
ously by Swift, a quadratic relationship provides a better 
fit than the linear fit, with of 81.8 and 162.0, respec¬ 
tively, for 11 DOF. The 0.3-3 keV flux and the > 200 


GeV flux are highly correlated, with a Pearson correla¬ 
tion coefficient (r) of 0.958. For the 3-7 keV band, the 
quadratic function fits the data better than the linear 
function, with of 58.0 and 114.0, respectively, for 11 
DOF. The r-value for the 3-7 keV flux as measured with 
Swift and the > 200 GeV flux is 0.954. 

6. MODELING THE BROADBAND SPECTRAL ENERGY 
DISTRIBUTION 

Previous MWL campaigns on Mrk 501 have been 
sufficiently characterized with a one-z one SSC model 
(|Acciari et al.ll201^ lAbdo et'ffi] 12011^ . although there 
are a few notable instances where a one-zone SSC model 
was found not to be appropriate for the br oadband emis¬ 
sion (jPian et al.ll998tTKataoka et al.ll999ll . In this study 
we decided to use the simplest approach, which is pro¬ 
vided by a leptonic model with a single emitting re¬ 
gion. The broadband spectral data were modeled with 
an equilibrium versi o n of th e sing le-zone SSC mo d el from 
iBottcher fc Chiaii^ (1200211 and iBottcher et al.l (I2013I1 . 
This model has been used to describe the broadband 
emission from various other VHE-detected blazars (e.g . 
Accipi et al.ll2009^ lAbdo et'ffi]l2011t^ lAliu et aLll201lL 

201311 . 

Within this equilibrium model, the emission originates 
from a spherical region of relativistic leptons with ra¬ 
dius R. This emission region moves down the jet with a 
Lorentz factor P. We set the Doppler factor J to 15 for 
all model representations. Notably, it has been shown 
that when using least-squares fitting of emission models 
to broadband data of Mrk 501, the Doppler facto r can 
vary widely from state to state (|Mankuzhivilll201^ . We 
do not complete least-squares fitting in this work and 
instead choose to fix the Doppler factor to 15 for the 
representation of all states, limiting the number of free 
parameters of the SSC model. The Doppler factor of 15 
is similar to the Doppler factor used in previous stud¬ 
ies of Mrk 501 (lAcciari et al.l 120111 : lAbdo et'all l2011al : 
iMankuzhivill 1201^ . In order to reduce the number of 
free parameters, the jet axis is aligned toward the line of 
sight with the critical angle 6 = 3.8°. At the critical an¬ 
gle, the jet Lorentz factor is equal to the Doppler factor 

Within this emission model, relativistic leptons are in¬ 
jected into this emission region continuously according 
to a power-law distribution Q = Qo'y~‘^ between 7 min 
and 7 max- The injected population of particles is allowed 
to cool. The simulation accounts for synchrotron emis¬ 
sion due to a tangled magnetic held Bq, Compton up- 
scattering of synchrotron photons, 77 absorption and 
the corres ponding pair production rates (via the general 
solution in iBottcher fc Schlickeiseiiri997[ 1. The cooling of 
the injected electrons is dominated by radiative losses, 
which are balanced by injection and particle escape from 
the system. This particle escape is characterized with 
an escape efficiency factor 7 = 100 , where Csc = rjR/c. 
The use of 7=100 is motivated by success in representing 
SEDs of Te V blazars in prev ious studies using the same 
model (e.g. lAliu et aT1l2013[) . The electron cooling rates 
and photon emissivity and opacity are calculated using 
similar ro utines of the code for j et radiation transfer de¬ 
scribed in iBottcher et al.l (|1997li . Together, the particle 
injection, cooling and escape mechanisms lead to an equi¬ 
librium particle population. 
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Figure 10. NuSTAR X-ray photon flux versus simultaneous > 200 GeV flux from MAGIC and VERITAS. The dotted lines show quadratic 
fits to the data, while the dashed lines show linear fits to the 3-7 keV and 7-30 keV bands. 
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Figure 11. Swift 0.3-3 and 3-7 keV X-ray photon flux values versus simultaneously measured > 200 GeV flux from MAGIC and VERITAS. 
The black and blue dotted lines show quadratic fits to the 3-7 keV and 0.3-3 keV data, respectively, while the black and blue dashed lines 
show linear fits to the 0.3-3 keV and 3-7 keV bands, respectively. For completeness, we also compare the linear and quadratic fits of the 
simultaneous 3-7 keV NuSTAR and > 200 GeV flux from MAGIC and VERITAS summarized in Figure 10 (light grey dashed and dotted 
line). 
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A key result of the equilibrium that occurs between 
continual particle injection, particle escape and radia¬ 
tive cooling is a break in the electron distrib ution yt, 
(referred to as 7 c within iBottcher et alJ 1201,111 , where 
tp. sc = trnniilh)- As des crjbed in Equations (1) and (2) 
of IBottcher et al.l (j 201 ,Hl . if 7 ;, is smaller than 7 ^™, the 
system will be in a fast cooling regime. If 7 ;, is greater 
than 7 mira, the system will be in a slow cooling regime. 
Within the fast cooling regime, the equilibrium particle 
distribution is a broken power law, with an index of 2 for 
particles with Lorentz factors less than 7 mm, and an in¬ 
dex of (g+l) for Lorentz factors above 7 mira- In the slow 
cooling regime, the resulting broken power law of the 
equilibrium particle distribution is equal to the injected 
spectrum [q) for particles with Lorentz factor below 7 ;,, 
and (g+l) above 7 b. It is known that a hard injected elec¬ 
tron spectrum would lead to a small amount of pile-up, 
followed by a smooth cut-off toward the high-energy end 
of t he distribution (for d etails, see, e. g. lKardesh^ll962l 
and lStawarz et all 1200811 . More specifically, the equilib¬ 
rium electron spectrum slightly deviates from the (9 + 1) 
approximation at the high-energy end (7 ^ 7 max) due to 
pile-up effects that increase as the injected spectrum be¬ 
comes harder (i.e. q < 1.5). Notably, although scattering 
in the KN regime is appropriately accounted for within 
the SSC model, neither the pile-up at the highest energy 
nor the energy loss (Compton cooling) of the electrons 
participating in scattering within the KN regime is ac¬ 
counted for within the model. The two aforementioned 
effects, however, are expected to result in a negligible 
deviation of the equilibrium electron spectrum from the 
approximated index of g -I- 1 . 

Lf. is the kinetic power in the relativistic electrons 
and As is the power in the Poynting flux carried by 
the magnetic field of the equilibrium particle distribu¬ 
tion. The Lf. and Lb parameters allow the calculation 
of the equipartition parameter Ls/Le- A state with 
an equipartition near unity minimizes the total (mag¬ 
netic field -I- particle) energy requirement to produce a 
given synchrotron flux. Therefore, from an energetics 
point of view a situation near equipartition is usually fa¬ 
vored. If the jet is powered by a Blandford-Znajek type 
mechanism, it is expected to be initially Poynting-flux 
dominated, and this luminosity is then (through an un¬ 
known mechanism, possibly magnetic reconnection) con¬ 
verted partially into particle energy. This conversion is 
expected to stop at an approximately equipartition situ¬ 
ation as an equilibrium is reached between the conversion 
of magnetic energy to particle energy, and vice versa (via 
turbulent charged-particle motion generating small-scale, 
turbulent magnetic fields). For ex amples of blaz a r mod ¬ 
eling based on equi partition, see iCerruti et al.l (j201,Hl : 
iDermer et al.l ()2014ll . Alternatively, a sub-equipartition 
magnetic field may be expected in an MHD-driven, ini¬ 
tially particle-dominated jet, where magnetic fields could 
be self-generated (amplified) by, e.g., shocks. The sub- 
equipartition magnetic fields that are often found in 
blazar SED modeling might therefore favor this latter 
scenario. Sub-equipartition states are a common result 
in the application of single -zone SSC emissio n scenar ¬ 
ios to VHE blazars as in lAliu et al.l (I^012alf3 12011 ); 
lAcciari et al.l (|2009allbl H l2008bll and lAbdo et al.l (|2011c l. 

The broadband data and model representations for five 
days from the MWL observation campaign are shown in 


Figure 12. The flux resulting from the model simula¬ 
tion (solid line) is corrected for absorption by interaction 
with the extragalactic background light (EBL) for the 
re dshift of z = 0.03 , assum ing the EBL model outlined 
in iDommguez et all (j2nilf l . The model thus represents 
the observed VHE emission as opposed to the intrinsic 
VHE emission. When applying the model to the data, 
the radio flux is likely to include a significant portion 
of extended radio emis sion and is therefor e taken as an 
upper limit, as done in lAbdo et al.l (j 2011 a[) . 

The parameters used to represent the data with the 
equilibrium model are summarized in Table 6 . The 
data in this work are represented with the emission 
model within the fast cooling regime, where the emit¬ 
ting equilibrium particle population follows n(e) oc 7 “^ 
for 7 t, < 7 < 'ymin and n(e) cx 7 “(‘?+i) for 7 mm < 7 < 
Imax- A particle population with an injection index of 
q = 1.8 — 1.9 provides a reasonable representation of 
the synchrotron emission on MJD 56395 (red; top panel) 
and 56420 (green; second panel from top). There are 
no Swift data for observations on MJD 56420. Each of 
these epochs (MJD 56395 and 56420) can be sufficiently 
described with similar parameters, although the SED on 
MJD 56420 requires a slightly more energetic electron 
population and lower magnetic field to account for the 
marginally elevated X-ray and VHE emission as com¬ 
pared to what is observed on MJD 56395. 

Although the highest VHE gamma-ray (> 200 GeV) 
flux during this campaign was observed by MAGIC on 
MJD 56484, a reliable spectrum from that MAGIC ob¬ 
servation could not be reconstructed due to the presence 
of Calima and the lack of LIDAR data to correct for it. 
Swift XRT also recorded the highest X-ray flux in its 
observation on the same day. On the other hand, there 
are sufficient broadband data to model the SED on MJD 
56485.0 (turquoise; middle panel Figure 12), which is less 
than one day later than the MAGIC and Swift observa¬ 
tion of the highest fluxes occurred. 

The light curve in Figure 5 shows that Mrk 501 dis¬ 
played relatively steady emission in each band between 
MJD 56420 and the elevated state observed by Swift 
and MAGIC on MJD 56484. In moving from the rel¬ 
atively quiescent SED on MJD 56420 to the elevated 
state observed on MJD 56485, a hardening of the in¬ 
jection spectrum is required (g=1.3) to match the X-ray 
spectrum observed by Swift XRT. With the injection in¬ 
dex responsible for the hardness of the synchrotron emis¬ 
sion at X-ray energies, the frequency at which the syn¬ 
chrotron emission peaks, is related to the spectrum of the 
injected particle population, and the magnetic field (Bq). 
When moving from the state on MJD 56420 to 56485.0, 
the strength of the magnetic field decreases, moving the 
peak of the synchrotron emission to lower energies. The 
decrease of the synchrotron flux resulting from a lower 
magnetic field is counteracted with an increase of particle 
luminosity Lg. Finally, to match the relative magnitudes 
of the synchrotron and inverse-Compton peak fluxes, the 
electron and photon density of the emission region was 
increased with a decrease of the emission region size. 
The decrease of the emission region size to 5.0 xlO^® 
cm on MJD 56485.0 provides a higher inverse-Gompton 
flux while maintaining the synchr o tron f lux. 

Following iBlumenthal fc Gouldl (jl97r)ll . the regime at 
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Figure 12. Observed broadband SEDs of Mrk501 on each of the days where NuSTAR observations occurred (red, green, blue and pink 
data). Additionally we include observations from MJD 56485.0 (turquoise, center panel), which show the SEE) one day after the most 
elevated flux state observed during this campaign. The broadband data are represented with a single-zone SSC model (solid line), with the 
model parameters summarized in Table 6. The Fermi LAT limits shown in the top two panels are taken from analysis of data between 
MJD 56381 and 56424, while the bottom three panels show Fermi results produced from analysis of data between MJD 56471 and 56499. 
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Table 6 

Single-zone SSC model parameter values (see Section 6 for overview of model and parameters). Model representations are shown along 

with the broadband data in Figure 12. 


Parameter 

MJD 56395 

MJD 56420 

MJD 56485.0 

MJD 56485.9 

MJD 56486.9 

7min [XIO^] 

2.0 

2.1 

2.0 

2.0 

2.0 

Imax [XIO®] 

1.0 

1.4 

1.4 

1.7 

1.4 

'Ibreak 10^] 

4.1 

4.6 

2.8 

3.3 

3.4 

9 

1.9 

1.8 

1.3 

1.3 

1.3 

V 

100 

100 

100 

100 

100 

5 

15 

15 

15 

15 

15 

Bo [G] 

0.06 

0.05 

0.03 

0.03 

0.03 

r 

15 

15 

15 

15 

15 

R [xlO^® cm ] 

7.0 

7.0 

5.0 

7.0 

7.0 

9 [degrees] 

3.8 

3.8 

3.8 

3.8 

3.8 

tvar [hr] 

4.3 

4.3 

3.1 

4.3 

4.3 

Le [xlO'*^ erg s“i] 

9 

12 

36 

28 

26 

e=Fs/Le 

1.8x10-2 

6.1x10-2 

5.3x10-4 

1.3x10-3 

1.4x10-3 
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which the up-scattering is occurring can be estimated 
(in the observer frame) according to 4 hi/gynpkT/^We c^, 
where 7 represents the energy of the electrons up- 
scattering z/synpk photons. If this quantity is less than 
1, the inverse-Compton emission is occurring within the 
Thomson regime, while if it is greater than 1, the emis¬ 
sion is occurring in the KN limit. With z^synpk at ap¬ 
proximately 5 keV, 4 ft.z/synpkbmin/'^We ^ 25, indicat¬ 
ing that, according to the model applied within this work, 
the inverse-Compton scattering of the photons near the 
synchrotron peak is far into the KN regime. We note 
that this is not necessarily in conflict with the quadratic 
relationship between the simultaneous X-ray and VHE 
flux measurements, but it implies a reasonably steady 
value of magnetic field which is supported by our SSC 
models; see Table 6 . For a more extensive discussion, see 
iKatarzvhskil (j2005[l . 

The SEDs on the days MJD 56485.9 and 56486.9 are 
similar to MJD 56485.0. All model representations ex¬ 
plored here result in emission scenarios which are heavily 
matter dominated (far below equipartition), where the 
majority of the energy is distributed within the particle 
population instead of in the magnetic field. Notably, even 
a single-zone SSC model is difficult to constrain, and the 
solutions presented here are not applied with the intent 
of constraining parameter space, but instead to just show 
that a reasonable representation of the data is possible. 
There are additional models (e.g. multi-zone or hadronic 
models) which might alternatively be used to describe the 
broa dband emission from Mrk 501 during these epochs 
(e.g. iTavecchio et al.l[2011l; lAleksic et al.ll2015bll . How¬ 
ever, these models have twice as many free parameters as 
single-zone leptonic models and, in this particular case, 
there are not strong constraints from MWL flux evolution 
correlations that point to the necessity of such models. 

7. DISCUSSION AND CONCLUSIONS 

The inclusion of the hard X-ray telescope NuSTAR in 
this observational campaign has provided unprecedented 
insight into the temporal evolution of the 3-30 keV X- 
rays emitted by Mrk 501. Before this campaign, Mrk 501 
had not been observed to display hard X-ray variability 
on timescales of ~7-hours. The fractional variability of 
Mrk 501 observed during this campaign was highly sig¬ 
nificant for the NuSTAR 7-30 keV band (Evar=0.6). 

Investigation of the DCF allows insight into possible 
leads or lags between the low (0.3-3 keV) and high (3-7 
keV) X-ray and VHE emission. The variability between 
these two bands shows evidence for a zero day lag. Cor¬ 
relation between the X-ray and VHE bands is further 
supported by the correlated variability inferred from the 
Pearson coefficients of 0.958 and 0.954 for simultaneous 
observations (occurring within one hour), respectively. 
Correlation is also found between the NuSTAR X-ray 
flux values and the simultaneous > 200 GeV flux val¬ 
ues (with observations occurring within one hour), with 
Pearson coefficients of 0.974 and 0.979 for the 3-7 keV 
and 7-30 keV bands, respectively. 

Correlation of variability between the X-ray and VHE 
flux, and more notably direct correlation without any 
lead or lag time, is a natural signature of SSC emis¬ 
sion. Within the single-zone SSC paradigm, the inverse- 
Compton flux is emerging from the same region as the 
synchrotron emission, and is fundamentally derived from 


the same particle and photon populations as the syn¬ 
chrotron emission. In this way, any variability in the syn¬ 
chrotron photon luminosity will immediately be trans¬ 
lated into a change in the up-scattered inverse Compton 
luminosity. 

In applying a single-zone equilibrium SSC model to the 
broadband data of Mrk 501, we find that the data could 
be reasonably represented in each of the five simultane¬ 
ous epochs. Notably, the injected particle populations on 
MJD 56485.0, 56485.9 and 56486.9 are very hard, with 
an injection index oi q = 1.3. Such a hard injection index 
is difficult to produce with standard shock acceleration 
scenarios alone, but is possible thr ough a magnetic recon¬ 
nection event (e.g. as explained i niRomanova fc L ovelacd 
Il992t iSironi fe Spitkovsl^ 12014 iGuo et al.l 1201411 . The 
increase in energy of the particle population (with an 
additional hardening to the injection index of q = 1.3) 
between the SED derived for MJD 56485.0 as compared 
to MJD 56420 indicates an introduction of additional 
energetic particles to the emission region, requiring some 
source of energy input. The decrease of the magnetic 
field, similar to what would naturally occur after a mag¬ 
netic reconnection event, is capable of accelerating par¬ 
ticles near the point of reconnection and producing the 
newly injected g = 1.3 particle population. Addition¬ 
ally, the decrease in the emission region size is consistent 
with a magnetic reconnection event that affects a more 
localized region as compared to a larger, more steady 
non-thermal emission region. More information on par- 
ti cle acceleration via rn agne tic reconnection c an be found 
in iWerner et al.l (1201411 and iGuo et all (1201511 . 

The variability timescale for these model representa¬ 
tions, quoted in Table 6 , is determined from the light¬ 
crossing timescale of the emission region according to 
Gar = R/c5(\ -I- z). For the emission region sizes and 
Doppler factor of <5=15 used within the model, the pre¬ 
dicted variability timescales of a couple of hours are 
compatible with the variability timescale observed dur¬ 
ing the broadband observations. The radiative cool¬ 
ing timescale is approximately equal to the synchrotron 
cooling timescale, tgync ~ 1-4 x 10"^ (Ho/0.06 G)“^7f7^ s, 
where 75 = 7 /( 10 ®). With a minimum light crossing 
time, corresponding to the minimum variability timescale 
of Gar 1.6 X 10“^ s (in the observer frame), all but the 
most energetic electrons within the emitting region cool 
on timescales that are longer than the crossing timescale, 
showing that the observed variability is likely a reflection 
of changes in the particle acceleration and/or injection 
processes directly. 

Notably, faster variability timesc ales have been ob- 
served from Mrk 501 in the past ('e.er. lAlbert et al.ll2007[) 
and so the model parameters shown here cannot be gen¬ 
eralized to all Mrk 501 flux variability episodes. NuSTAR 
observations show the hard X-ray flux to significantly 
decrease by more than 10% between its 90-min orbits. 
Moreover, on MJD 56420 the source hard X-ray flux was 
observed to change by a factor of greater than 40% in 
the 7-30 keV band during a 7 hour exposure. 

In an attempt to describe a possible emission scenario 
which might result in the broadband SED variability ob¬ 
served for Mrk 501 in 2013, the parameter changes were 
made to the single zone equilibrium SSC model monoton- 
ically. With a degeneracy between several of the input 
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parameters, the model applied here cannot be used for 
conclusive studies regarding which changes occur within 
the emitting region from one state to the next. Instead, 
through the study of band-to-band spectral variability, 
leads and/or lags and fractional variability, as well as 
broadband modeling of various flaring episodes, we find 
compelling evidence to support a single zone SSC emis¬ 
sion scenario for Mrk 501 during the broadband obser¬ 
vations in this campaign. 

The collection of simultaneous broadband observations 
is a necessity for the study of the relativistic emission 
mechanisms at work within blazars such as Mrk 501. It 
is known that these sources vary continually, with charac¬ 
teristics that significantly change between different flar¬ 
ing episodes, requiring the continuation of deep broad¬ 
band observations such as those presented in this work. 
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